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Abstract 

The Green’s function (GF) for the steady state Laplace/Poisson equation is derived for an 
anisotropic finite two-dimensional (2D) composite material by solving a combined Boussinesq- 
Mindlin problem. The source term for the GF is a delta-function located somewhere in the bulk 
of the solid (Mindlin problem). The boundary condition for the GF is prescribed in terms of a 
delta function at a single point on the outer boundary (Boussinesq problem). The calculated GF 
is, therefore, quite general and can be used to find the full solution of the Faplace/Poisson 
equation for an arbitrary but integrable distribution of sources and boundary values in a 
composite material. The piecewise continuous composite, consisting of a host material 
containing one or more inclusions, is assumed to be infinite and periodic in one direction (X- 
direction) but of finite arbitrary size in the perpendicular direction in the plane (Z-direction). 

A semi-discrete semi-analytic model is developed for the calculation of the GF subject to the 
boundary conditions described above. The primary equation is solved by using a partial Fourier 
transfonn technique. Our model is semi-discrete in the sense that the solution space is discretized 
only in the periodic X-direction but not in the Z-direction. This allows the space to be filled up 
exactly in terms of the discrete elements even if the shape of the inclusions are geometrically 
irregular. Further, our model is semi-analytic in the sense that the Fourier integration in the Z 
direction is done exactly. Hence a 2D problem needs only a ID discretization. An approximate 
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analytical estimate shows that the numerical convergence of our semi-discrete model is at least 
an order of magnitude faster than conventional fully discretized models (e.g., finite-element or 
boundary-element models). 

As an application of the derived method, GF numerical results are presented at all points on the 
plane of a 2D nanocomposite solid, consisting of a phosphorene monolayer containing an array 
of metallic inclusions. This method will be useful in modeling heat flow and electrostatic 
potential distribution in modern 2D anisotropic material systems. These problems in general, and 
phosphorene-based material systems in particular, are of strong topical interest in view of their 
possible application in solid state devices for energy conversion and quantum computing. This 
paper is another step towards developing characterization techniques of modern 2D materials 
using the GF, as suggested in an earlier paper [1]. 


1. Introduction 

Modern two-dimensional (2D) materials are being considered as primary materials for 
thennoelectric and quantum computing devices and many other revolutionary applications [2] . 
Currently there is a special interest in phosphorene because of its unusual physical properties and 
possible application in solid state devices [3-8]. Its thennal conductivity is highly anisotropic, 
which makes it especially interesting, as well as relatively more challenging, for mathematical 
modeling. For example, the thennal conductivity of black phosphorene has been estimated [9] to 
be 110 W/m-K and 36 W/m-K along its armchair and zigzag directions at room temperatures. 
Modeling and characterization of these materials are, therefore, subjects of strong topical 
interest. 

An earlier paper [1] suggested that the Green’s function (GF) is not just a mathematical artifact 
but is a physical quantity that can be measured using a point probe technique such as SPM 
(Scanning Probe Microscopy). Mathematically, GFs are known to be useful tools for solving the 
Laplace/Poisson equation. If they can be measured, they can provide a direct and reliable 
characterization technique that will better enable the industrial applications of modern 2D solids. 
However, to deconvolve the measurements and to extract physically meaningful data, it is 
necessary to have highly efficient and robust computational tools for calculating the GFs and for 


2 



identifying the discriminants in the GF that provide the characteristic features of the materials. 
This paper is a step in that direction, with special reference to phosphorene and its composites. 

Formally, the GF [10, 11] gives the response of a system to a point probe, which is exactly what 
is measured by SPM type techniques. At macroscales, the 2D GF for an infinite solid in the 
continuum model can be written in tenns of logarithmic functions. For finite solids, the GF must 
be calculated by specifying appropriate boundary conditions. Such problems can be broadly 
classified as the classic Boussinesq or Mindlin problems [11-13], In the Boussinesq problem, 
Dirichlet boundary conditions are prescribed at the boundary: the value of the unknown function 
is forced to unity at a single point on the boundary and zero at all other points on the boundary. 

In the Mindlin problem the point probe, referred to as the source, is buried in the bulk. The 
Mindlin and the Boussinesq GFs can then give the solution of problems of physical interest for 
an arbitrary distribution of sources subject to prescribed boundary conditions using linear 
superposition. For brevity, we will refer to this problem as the BM (Boussinesq/Mindlin) 
problem and the coupled GF for this problem as the BMGF. 

A BM problem requires solution of an inhomogeneous equation (Poisson equation) subject to 
inhomogeneous Dirichlet boundary conditions. Our approach separates out the inhomogeneity in 
the boundary condition. We then need to solve an inhomogeneous equation subject to 
homogeneous boundary conditions, which gives the primary solution. The boundary values are 
reproduced by a separate integral over the boundary. We use the complex Cauchy representation 
of the singular GF integral. The real part of the integral gives the primary solution whereas its 
imaginary part is used for generating of the boundary values. Integration over the boundary is a 
common procedure in the conventional GF method [10]. However, a relative advantage of our 
technique over the conventional method is that we don’t need to evaluate the derivative of the 
GF. Calculation of the derivative of the GF is not always straightforward because of the 
discontinuities in the GF. 

In this paper, we calculate the BMGF for a general case then apply it to an anisotropic 2D 
bimaterial composite solid, with phosphorene as the matrix. We start with a model solid 
containing a unit point source in the bulk with boundary value of the unknown function 
prescribed to be unity at a single point on the boundary and zero everywhere else. Both the host 
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and inclusion are assumed to be homogeneous and piecewise continuous. Our formulation is 
valid for the electrostatic and thermal response of any 2D material. 

The primary motivation for this work is to develop efficient techniques for characterizing 2D 
composite materials, for which measurement and modeling of electrical or thermal conductivity 
is an established technique [14-18]. Many papers are available in the literature on the solution of 
the Laplace/Poisson equation for composite materials, with excellent reviews [18-22] . In real 
cases, an analytical closed form solution is seldom available, so that it becomes necessary to 
resort to numerical techniques for solving the equations. In such techniques, problems of 
stability, convergence, and possible errors associated with the approximations must be 
considered. Accelerated convergence of thermal modeling is crucial in problems of materials 
design [23]. The possibility of measurement of our GF, as suggested in our earlier paper [1], is a 
step towards achieving this objective. 

In many cases of physical interest, the inclusions can be modeled as fonning a regular array in 
the solid. These cases can be simulated by specifying periodic boundary conditions in an infinite 
solid [24]. For such systems, the Fourier transfonn provides a particularly useful representation 
of the GF [10-12, 25-27]. Since the GF is essentially an operator, its calculation requires a 
representation in terms of mathematical functions [28]. The Fourier representation is a powerful 
representation of the GF because it can be used for representing many different mathematical 
functions of physical interest. The computational cost of the GF technique arises largely from the 
integration for calculating the inverse Fourier transform. The calculation of GF for a 2D material 
will require a 2D integration of a singular integrand [28], In general, the Fourier integrals need to 
be evaluated numerically except for the simplest material geometries. In this paper, we describe 
a partial Fourier transfonn (PFT) technique for calculation of the GF. This technique is 
convenient for analytical simulation of the periodicity of the solid in specific directions. 

We have developed a semi-discrete model for calculating the GF for 2D Laplace/Poisson 
equations for thennal or electrostatic simulations. Our method accounts for anisotropic 
conductivity. Although many papers have been published on modeling the thermal/clcctric 
response of 2D nanomaterials, very few papers are available in the literature on anisotropic 
thermal/electric GFs. Elegant integral equations for anisotropic GFs have been derived by 
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Chang [29] for thermal applications and by Hanson [30] for interaction of electromagnetic waves 
with 2D materials, but were not applied to 2D composite materials of interest to us. 

The GF method is an alternative to the more common numerical techniques for solving boundary 
value problems, which are based upon the use of purely numerical finite and boundary element 
methods (FEM, BEM) and can be used to solve elliptic equations [14-16, 31-36]. The primary 
computational advantage of our GF technique, presented here, is that it is partly analytical. Even 
in the general case, the integration over one of the Fourier variables can be done analytically (see 
Appendix A), which results in a considerable economy of computational effort. This technique 
has the additional advantage that it can be seamlessly linked to atomistic models like the Bom- 
von Karman model and lattice statics GF [28, 37], which is useful for multiscale modeling. A 
somewhat subtle but important advantage of the GF technique for nanomaterials is that the 
limiting cases like infinite solid or point charges in a real solid can be simulated precisely. This 
would avoid spurious size effects, particularly in dynamic problems, which are inherent in purely 
numerical techniques. 

This paper is the first in a series of papers on the formulation and application of the semi-discrete 
GFs for composites. In this first paper, we will derive the BMGF for a 2D nanocomposite 
consisting of a metallic inclusion in phosphorene. These inclusions can be used to tune the 
effective conductivity of the composites. The effect of inclusions on thennal and electrostatic 
properties has been extensively studied in graphene [38] but not in phosphorene. In this paper, 
we assume that the metallic inclusion has infinite conductivity relative to that of the host solid. 
This serves as a limiting case for general inclusions. Later papers in this series will deal with 
inclusions of zero and intennediate conductivities and measurable average conductivities of the 
composites. 

In this paper, our interest is confined to monolayer or atomistically thin 2D material systems. An 
important assumption inherent in the use of Poisson/Laplace equations for thermal/elastostatic 
modeling is that the continuum model is valid for the material and discrete atomistic effects can 
be neglected. This assumption is common and should be valid at length scales much larger than 
the interatomic spacing in the 2D solid. Note that phosphorene and many other 2D solids are not 
exactly planar and do have a finite width in the third dimension. However, they can be treated as 
perfect 2D planar solids at the length scales when Poisson/Laplace equations can be used. Note 
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also that the ‘discreteness’ of our semi-discrete model is a computational assumption and has 
nothing to do with the atomistic discreteness of solids. 

In the literature, the tenn Green’s function has been used in different meanings in different 
contexts. It can be an operator, which is defined by its effect on mathematical functions or it can 
be a mathematical function itself, which is a solution of an operator equation. The main 
difference between the two is that an operator does not need boundary conditions, whereas a 
solution function does. In the latter case, it is also called the response function. Primarily, the GF 
is the response of a solid to a single point source and can be represented by an operator as well as 
a mathematical function. The responses are additive and the total response of the solid to a 
distribution of point sources is obtained by adding the responses of the individual sources. 

In our present case, we have two different kinds of point sources that are non-additive. One 
source is buried in the solid, which we call the Mindlin source. The other source arises from the 
boundary value prescribed at a single point on the outer boundary, which is the Boussinesq 
source. In this paper we calculate the response function of the solid to these two sources subject 
to all the prescribed boundary conditions. The two sources are coupled and the net response of 
the solid is not a linear sum of the responses to the two sources. Of course, sources of the same 
kind are additive among themselves. For example, the responses of a distribution of Mindlin 
sources for the same Boussinesq source can be added together. Same applies to a distribution of 
Boussinesq sources. The advantage of our approach is that the two sources can be varied 
independently. This makes our model more efficient for simulating a much wider class of 
material systems. This will also make the measured GF more useful for materials 
characterization. 

In Sec. 2, we describe the main features of our semi-discrete model and the use of the PFT 
technique for calculating the BMGF of a homogeneous solid. The GF is calculated in modules 
that we build up in successive stages of complexity. The homogeneous solid is our basic module, 
and we introduce inclusions in Sec. 3. The application of the method to phosphorene containing 
a metallic inclusion is described in Sec. 4. Finally, a summary of conclusions is presented in Sec. 
5, along with general discussion. 
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2. BMGF for a homogeneous solid 

The model solid geometry is shown schematically in Fig 1 along with a Cartesian frame of 
reference. The figure actually shows a composite but for the purpose of this section, we consider 
only the host region labelled as H and assume that the inclusion I is identical to H. The 2D plane 
of the solid is assumed to be the XZ plane. The X and Z coordinates will be denoted by x and z 
respectively and a vector on the XZ plane will be denoted by r = (x, z), with r t = x and r 3 = z. 

The Cartesian components will be labelled by the indices i, j (=1 or 3), with no reference to the 
Y-coordinate in this 2D paper. Summation over repeated Roman indices will be assumed unless 
stated otherwise. 

The host material is assumed to be a rectangular continuum. It extends from z= -d to z=+d in the 
Z-direction and from x=-L to x=+L in the X-direction with 2L as the period in the X-direction. 
We neglect the discrete atomistic structure of the material, so our calculations are valid at length 
scales larger than the atomistic dimensions. This is consistent with the continuum approximation 
inherent in the Laplace/Poisson equation. It allows us to neglect the atomistic zig-zag structure of 
the edges and any ripples or unevenness on the plane. Thus, the treatment given here is 
applicable to general 2D materials. 

In an experimental set up, a probe such as an electrode or a source of heat can be placed on the 
solid. The resulting potential or temperature distribution can then be measured across the surface 
of the solid. For a 2D solid, in principle, this gives the GF [1], We simulate this system by 
assuming cyclic boundary conditions in the X-direction and Dirichlet boundary conditions in the 
Z-direction. In the X-direction, our model in Fig. 1 can be identified as a super-cell of the Born- 
von Kannan model [37] of a crystal lattice. 
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Fig. 1: A two-component composite solid in 2D. The solid is confined to the Cartesian XZ 
plane as shown in the figure. The solid consists of a metallic inclusion PQST, labelled I, 
inside the host solid ABCD, labelled H. The source point, where the probe is applied is 
denoted by M and the field point, where the response is calculated/measured is denoted by F 
Both F and M are assumed to be in the region H, which includes the boundaries ABCD and 
the interface PQST. 
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In the GF method, we formulate the problem in modules of increasing complexity by imposing 
appropriate boundary conditions. In the BEM literature, the zeroth order module is usually an 
infinite homogeneous solid and the GF is referred to as the fundamental solution. Our zeroth 
order module is also an infinite homogeneous solid, but it is periodic in the X-direction. To 
prescribe the cyclic boundary conditions in the X-direction, we assume a semi-discrete model of 
the solid. In this model, the X- component of r in the solid is allowed only discrete values over 
the range -L to L, whereas the z-component of r is assumed to be a continuous variable that 
extends from - oo to + oo. In the next module, we will restrict z to the range -d < z < d by 
prescribing boundary conditions at z = -d and d. 

The temperature or the potential <D(r) in an anisotropic medium is the solution of the Poisson 
equation, which is our primary equation 

X O(r) = - p(r), (1) 

where the operator X is defined as 

X=V.CV, (2) 

V= (d/dx,d/dz) (3) 

is the 2D gradient operator on the XZ plane, C is the second rank conductivity (electrical or 
thennal) tensor of the material, and p (r) is the source density. For electrostatic problems, a 
point charge is the source, whereas in thermal problems the source can be a source of localized 
heating. In what follows, for brevity, we will refer to O as the solution function. 
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In terms of the Cartesian indices, Eq. (2) can be expressed as 

X — VfCyVy. (4) 

For an isotropic solid, the conductivity is scalar, and the tensor C is of the form C q. In this 
case, the operator x reduces to the conventional form: 

X — C V 2 (for isotropic solids). (5) 

We prescribe the following boundary conditions on <D(r): 


i. Cyclic boundary conditions on the X-axis: 

0 (x+2Ln, z) = O (x, z), (6) 

where 2L is the length of the supercell or the periodicity in the X direction (see Fig. 1), and n is 
any integer. In view of Eq. (6), we can assume that the solid is infinite in the X-direction and x is 
confined to the unit cell, that is. 

-E < x , x < E. (7) 

ii. Function at the top and the bottom edges: 

0(x,z = d) = V d (x), (8) 

0(x, z = —d) = 0. (9) 

where E d (x) is a kn own function. For example, it can be the applied potential at the top surface 
(z=d) in an electrostatic problem or the prescribed temperature in a thermal problem. Eq. (8) is 
an inhomogeneous Dirichlet boundary condition for the Boussinesq problem [11]. 

The above equations define the classic boundary-value problem. We need to solve the 
inhomogeneous primary equation Eq. (1), subject to boundary condition given by Eq. (6) in the 
X-direction and the inhomogeneous boundary conditions given by Eq. (7) and (9). Equation (9) 
as such is of course homogeneous but our approach is valid for general boundary values, which 
may be different for different segments of the boundary. 
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It is convenient to express the solution for the inhomogeneous boundary condition in Eq. (8) as a 
sum of two parts as described below. We separate out the inhomogeneity in the boundary 
condition by writing the solution function as 


O(r) = O 0 (r)+ 0 B (r), (10) 


where <3>o is the solution of Eq. (1), corresponding to the homogeneous part of the boundary 
value, that is, Eq. (8) with its RHS=0, and Ob is a suitably chosen function that accounts for the 
inhomogeneity in the prescribed boundary value. We will refer to Oq as the primary solution and 
Ob as the boundary solution. We need to choose a function for the boundary solution, which has 
the following properties: 

O b (x,z = d) = V d (x) (11) 

and 

O s (x,z = -d) = 0. (12) 


Using Eq. (10) in Eq. (1), we obtain 


xO 0 (r) = -p eff(r), (13) 

where 

Pe//0) = P (r) + X^bO). (14) 

So far, the choice of O s (r) is arbitrary except that it satisfies Eqs. (13) and (14) We choose 
0 B (r) such that it is a solution of the homogeneous part of Eq. (1), which is 

X$ B (r) = 0. (15) 
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In this case 


Pe//0) = P (r). (16) 

Hence we need to solve Eq. (1) for homogeneous boundary conditions which would give <3>o i n 
Eq. (10). In many practical cases, the boundary conditions can be mixed type for different 
segments at the boundary, which can also be simulated with our GF approach. 

In the GF method, instead of solving Eq. (1) directly, we solve an equivalent operator equation 
for a delta function type inhomogeneity on the RHS. Further, we replace the inhomogeneous 
boundary condition in Eq. (8) by a delta function type inhomogeneity. A delta function 
inhomogeneity is located at a single geometrical point. It acts as a source point because it makes 
a well-defined and isolated contribution to the solution. The geometrical point at which the 
solution is calculated is called the field point. The solution at the field point is referred to as the 
response of the medium to the stimulus applied at the source point. Thus the GF is a function of 
two points- the field point r, and the source point r , which will be the arguments of the GF. 

To avoid any confusion, we will separate the source and the field points in the argument of the 
GF by a semicolon, whereas the components of a vector will be separated by commas. The field 
point will be the first argument and the source point will be the second argument of the GF. Thus 
the GF at the field point r due to a source at r is denoted by G(r ; r ) or G(x,z ; x ,z ), 
where r and r are 2D radius vectors with components (x,z) and (x ,z ), respectively. In 
homogeneous solids which have translation symmetry, G(r ; r ) will depend upon r and r 
only through their difference and not upon them individually. In such cases, it is conventional to 
express G as a function of single argument r-r . Accordingly, if GF is expressed as a function 
of a single argument such as G(r), it implies that r =0. 

The mathematical problem that we need to solve for the GF is now defined as follows: 

Primary equation: 

xG(r; r' ) = - 6(r-r'), (17) 

where G (r; r ) obeys homogeneous boundary conditions, and S (r) is the Dirac delta 
function. For a scalar variable the delta function is defined by the following relations 

5(£) = 0 for f * 0 (18) 
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and 


/ 5(f) dr = 1. 


(19) 


where the integration is over all space. The delta function of a vector is defined as the product of 
the delta functions of its components: 


S(r ) = 5(?i) S(r 3 ). 


( 20 ) 


Note that X operates on r but not on r in Eq. (17). 

To account for the inhomogeneity at the boundary, we introduce a new function g B (r, r ), 
which is a solution of the following boundary equation 

Xg B (r;r') = 0. (21) 

Boundary conditions 

i. Periodic boundary conditions on the X-axis: 

G(x+2Ln, z; x , z ) = G(x, z; x , z ), 
g B (x+2Ln, z;x ,z )= g B (x, z; x ,z ), 

Boundary values at the top and the bottom edges 
G(x, z = ±d;x ,z ) = 0, 
g B (x,z = d;x ,z ) = 8(x — x ), 


( 22 ) 


(23) 


u. 


(24) 


and 


(25) 


(26) 


g B (Xj z = -d;x , z ) = 0, 

for all values of z (actually, independent of z ). In the above equations and in what follows, x 
and X are assumed to be confined to the unit cell as in Eq. (7). 

We can now write the primary and the boundary solutions as follows 

<J> 0 (r) = / G(r, r') p(r')dr', (27) 

where the integration is over all space, and 


(r) = g B (x,z; x' ; z') V d (x')dx', (28) 
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where the integration is over the x-values at the boundary. The total solution of Eq. (1) is then 
given by Eq. (10). It can be verified by direct substitution that the total solution satisfies Eq. (1) 
in view of Eqs. (17) and (21). It would will also satisfy the boundary conditions because G and 
g B satisfy Eqs. (22) - (26). 

Effectively we have expressed the solution of an inhomogeneous equation (Eq. (1)) subject to 
inhomogeneous boundary condition given by Eqs, (8) in two parts. One part is O 0 (^) , which is 
solution of an inhomogeneous equation subject to homogeneous boundary conditions and the 
other part is d> B (r), which is solution of a homogeneous equation subject to inhomogeneous 
boundary conditions. The boundary value is generated by integration over the boundary. This is 
somewhat similar to the boundary element method [10, 32], except that, in the present method 
the derivative of the GF is not needed for this calculation. 

The main computational effort is the construction of the G (r; r ) that satisfies the 
homogeneous boundary conditions. We write 

G(r;r )= G p (r;r ) + G h (r;r ). (29) 

where G p (r; r ) is the particular solution of Eq. (17) and G h (r; r ) is a solution of the 
homogeneous part of the same equation, that is, for the RHS=0. Thus the two terms on the RHS 
of Eq. (29) are defined by 

X G p (r; r' ) = - 8(r — r'), (30) 

and 

xG„(r;r ) = 0. (31) 

The homogeneous solution G h (r; r ) may not explicitly depend upon r . However, it may 
indirectly depend upon r because it contains arbitrary constants that are detennined by 
imposing the boundary conditions. The form of the functions G h (r; r ) may or may not be the 
same as that of g B (r; r'). Both are solutions of the same homogeneous equation. The difference 
is that the arbitrary constants in G h (r; r ) are determined from the homogeneous boundary 
conditions imposed on G(r; r ) whereas in g B (r; r') they are determined from the 
inhomogeneity in the boundary condition. 
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For the BMGF, we calculate the response of the solid to a point source of strength F M located at 
r M . = (x M ,, z M )■ We will refer to it as the Mindlin source. It amounts to solving Eq. (1) for the 
following source distribution 

p (r) = F m S(r - r M ). (32) 

In this paper we calculate the solution function for p (r) given in Eq. (32). The GF is defined as 
the response to a unit source. Accordingly, we assume F M = 1 in Eq. (32). The GF will be a 
function of r and r M but for notational brevity, we will not show its dependence on r M , except 
where it is needed for clarity. Using the distribution given in Eq. (32), we see from Eq. (27) 

O 0 ( r ) = G(r;r M ). (33) 

Note that G(r; r M ) and d> 0 (r) satisfy homogeneous boundary conditions at the outer 
boundaries as prescribed by Eq. (24). 

The boundary condition on the BMGF is prescribed to be as a delta function localized at a single 
point on the outer boundary. The boundary conditions are given by Eqs. (11) and (12), where 


V d (x) = V 0 8(x - x„), (34) 

The point (x B , d) acts as a second source of the GF, which is the Boussinesq source. 

The PFT representation of the GF is essentially an adaptation of the conventional Fourier 
representation. The general GF in the full 2D Fourier representation [1, 10, 11, 28, 32, 33, 39, 
40] is written as given below: 

G(r)= (1/2t r) 2 J_ + ” G(k) exp(tk. r) dk, (35) 


where l = V— 1 , k is a 2D reciprocal vector conjugate to r, and G(k) is the inverse Fourier 
transform of G(r). It is given by 


G(k)= (1/2 n) 2 J_ + “G(r) exp(-tk.r) dr. (36) 
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For free infinite space, 


—oo < kj , Tj < oo . (37) 

In deriving Eq. (36), we have used the orthogonality relations valid for each coordinate i 
separately (no sum over i): 

(1/2tt) C exp(-t) dr z = S(fc f )- (38) 


As described earlier, Eq. (35) implies that the source is at r =0. For a source at r M , we write 
using Eqs. (4) and (35), 

X G(r; r M ) = —(1/27 1 ) 2 J_ + “ <P(k) G( k) exp[tk. (r - r M \ dk , (39) 

where, for symmetric C, 

^(k) = (C n + 2C 13 k 1 k 3 + C 33 k\ ). (40) 

Using Eq. (39) in Eq.(30), and the orthononnality of the exponential functions, we obtain 

G( k) - H'(k)- 1 . (41) 

Eq. (41) shows that the integrand in Eq. (35) has pole(s) where *?(&) = 0. The integral is 
defined in the Cauchy sense and has been evaluated in Appendix A. 

The Appendix also describes the construction of a homogeneous GF using the same 
representation as for evaluating the Cauchy integral. For future reference, we note here that 
another independent homogeneous GF can be derived by using Eq. (30). This equation shows 
that if r is outside the solution space, the RHS is 0 for all values of r. This is because the 
argument of the delta function on the RHS cannot be 0 for r inside the solution space. This is 
the basis of the method of virtual sources or charges [1,41]. We will use it for constructing 
BMGF in the next section. 

Now we impose the cyclic boundary condition given by Eq. (6) in the x-direction. For this 
purpose, we treat x and ki as discrete variables with the following constraints: 

Range: - L < x < + L and (-7t/L) < k, < (n/L) (42) 
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Allowed values: ri — r ln and kj — k ln , where 

r ln = n L/N and ki n = n n/L, (43) 

N is a positive integer, and n is a positive or negative integer, including 0. In view of the 
constraint (42), 

-N < n< N. (44) 

Henceforth, for brevity, we will use the notation x for r in , z for and q for ki n without the 
subscript n. We will denote vectors in real and reciprocal spaces in terms of their components as 
follows: 

r = (x, z) and k = (q, k 3 ). (45) 

Thus, in our semi-discrete model, z and k 3 are ID continuous variables, whereas x and q are ID 
discrete variables in the vector space defined by Eq. (43). We note that, with these choices, our 
model becomes analogous to the Born - von Karman model [28, 37, 42] of a monoatomic 
crystal lattice with cyclic boundary conditions in the x-direction. The range of the q values in Eq. 
(42) can be identified with the first Brillouin zone of the lattice. Our formulation can then be 
recognized as the super-cell approach for modeling of crystal lattices with defects containing 2N 
atoms in a single supercell. Without the periodicity, N is formally infinite and our model reduces 
to infinite Fourier series solution of the Laplace/Poisson equation. These solutions have been 
discussed in detail in the classic literature [11, 22, 25, 26]. 

In this model, the integration over ki has to be replaced by a sum over the allowed values of q as 
follows: 


(1/2tt) f” f( k{) d fc, = (1/2N) 2, f(q) (46) 

for any function f. The corresponding orthononnality relations for the Fourier functions are given 
by: 


(1/2A0 X q exp (tqx) = 6 (x), (47) 
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and 


(1/2 N E x exp (iqx) = 5 (q). (48) 

The summations in Eqs. (47) and (48) are over the allowed discrete values of q and x given by 
Eq. (43). The delta functions in the above equations are functions of discrete variables as used 
in Born-von Kannan lattice models [37]. In numerical evaluation of the discrete lattice sums 
over q as well as x, it is important to account for the weight of each tenn. In a supercell, the 
weight is 1/2 for the end terms and unity for other terms. 

Using Eq. (46), we write Eq. (113) of Appendix A as follows: 

r; r M ) = 1 / 2N 'Z q exp[tq(x - x M )] T p {q,z - z M ) , (49) 

where, for any z, 

r p(q,z)= (^) p /_ + r [ ( ?(q^ 3 )]~ 1 exp(L/c 3 z)d/c 3 (50) 


Evaluation of the integral in Eq. (50) is given in Appendix A. The result, for general z, is 


T p (q,z)= exp(ip R qz - ^q\z\) /(2C 33 q^), (51) 

where /? has been defined in Appendix A. 

It can be verified directly that G p (r; r M ) given by Eq. (49) satisfies Eq.(30). It is, therefore, the 
desired particular solution in Eq.(27). Further, because of the choice of q-values, the solution 
satisfies the cyclic condition at x = 2L in Eq.(6) . Note that G p (r; r M ) depends upon r and r M 
only through their difference because the particular solution is for a perfect infinite solid that has 
full translation symmetry. 

To obtain the homogeneous solution, we note from Eq. (127) of Appendix A that only those 
values of (q, k 3 ) are allowed that satisfy ^(k) = 0. These are the same roots for k 3 that 
contribute to the particular solution. The only difference is that all the roots contribute to the 
homogeneous solution irrespective of their location in the complex k 3 plane, which can be in the 
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UHP (upper half plane) or LHP (lower half plane). The result, therefore, would not depend upon 
the sign of z. Thus, we obtain, after some algebraic manipulation: 


Gfc(r) = 1 / 2 n Zq exp(tqx) Y h (q, z) , (52) 


where 

T h (q,z)= [6>i(q) exp(t^z) +0 2 (q) exp(t^z)], (53) 

9i and 02 are arbitrary constants, which can be functions of q. These constants are detennined by 
imposing the boundary conditions at the top and the bottom surfaces of the solid given by Eqs. 

(8) and (9), respectively. We then have 

G( r; r M ) = G p (r; r M ) + G h (r), (54) 

where G p and G h are, respectively, given by Eqs. (49) and (52). As discussed in Appendix A, 

G h (r) does not explicitly depend upon r M . It will of course have an implicit dependence on r M 
through the values of the arbitrary constants after the two solutions are coupled through the 
boundary conditions as shown below. 

In terms of its PFT, we write the GF on the LHS (left-hand side) in Eq. (54) formally as 

G(r; r M ) = 1 / 2N exp(tqx) T t (q, z; z M ) . (55) 

From Eqs. (49), (52), and (54), we can write 

r tiq.z) z M ) = T p (q,z - z M )exp(-Lqx M ) + r h (q,z) (56) 


Using Eq. (47) and Eqs. (49) - (53) in Eq. (54), and the boundary conditions Eqs. (8) and (9), we 
obtain the following linear equations that uniquely determine the arbitrary constants 0i and 0 t: 

0i(q) exp(t/?qd) + 0 2 (q) exp(t/?qd) = -T p (q, d - z M ) exp (~iqx M ), (57) 
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and 


0i(q) exp(-i^qd) + 0 2 (q)exp(—i^qd) = 

-T p (q,-d-z M )exp(-iqx M ). (58) 

In writing Eqs. (57) and (58), we have used the fact that d > z M . 

In general, C is a symmetric tensor. It can be diagonalized by rotating the coordinate axes. In this 
representation, which is obviously different for different solids, Cn = C 31 =0. In this case, as 
shown in Appendix A, we have 


P = + Cl, 

(59) 


Pr = 0 and Pi = C, 

(60) 


where 



C 2 — Ci 1 /C 33 . 


(61) 

Using Eq. (126) (Appendix A), we have 



T p (q,z - z M ) = exp(— a? \z-z M \) 

/(2q V Cn 

£ 33 )- 


Instead of the exponentials, it is possible to choose hyperbolic sinh and cosh functions in Eq. 
(53), which may be computationally more convenient for finite values of d. We write r\ in 
terms of hyperbolic functions as follows: 

r ft (q, z) = [D x (q) cosh(cqz)/cosh (cqd)+ 

D 2 (q) sinh (cqz)/ sinh (cqd)], (63) 

where D 12 (q) are arbitrary constants, which like 0 1.2 in Eq. (53), are functions of q and will also 
depend upon r M . 
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Following the steps leading to Eqs. (57) and (58), we obtain the following relations for Di(q) and 


D 2 (q): 


M q) = exp (—t qx M ) r sM (q) , 

(64) 

D 2 (q) = exp(—t qx M ) r mM (q) , 

(65) 


r 5M (q) = [r p (<7, d - z M ) + r p (q, -d - z M )\/2 , (66) 

and 

^mM Cq) [^p (jj>d z M ) Fp (c[, d z m )\/2 . (67) 

Equations (64) and (65) have been derived in a frame of reference in which the conductivity 
tensor is diagonal, that is, C 13 = C 31 = 0. The above equations will of course reduce to the 
isotropic case for c = 1. In writing Eqs. (64) and (65), we have used Eq. (62) and the fact d > |z M | 
for a Mindlin problem. 

It can be verified by direct substitution and after a rather lengthy algebraic calculation that the 
expression given in Eq. (54) satisfies Eq. (17) for r = r M and the boundary conditions given by 
Eq. (24). In Eq. (54), T p and Th are, respectively, given by Eqs. (62) and (63) using Eqs. (64) 
and (65) for Di and D 2 . Eq. (54) is, therefore, the desired GF for homogeneous boundary 
conditions. 

By definition, GF is the solution of Eq. (1) for a unit source at r M or Fm = 1 in Eq. (32). For 
F m ^ 1, the particular solution gets multiplied by F M - Note that the boundary conditions will be 
satisfied for any value of F M or r M . 

Finally, to get the desired solution of Eq. (1) that satisfies the boundary conditions given by Eq. 
(11) and (12), we need to construct the boundary function g B (r; r ), which satisfies Eqs. (21), 
(23), (25), and (26). As shown below, we can construct g B (r; r ) from the homogeneous 
solution derived in Appendix A. 

First, to ensure the periodic boundary conditions given by Eq. (23), we simply use the same PFT 
representation as used for the primary GF G. We obtain from Eqs. (23) and (25): 
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g B 0,z;x > z )= 1 / 2N Iq exp[Lq(x-x )] g h (q,z) 


( 68 ) 


where 

g h (q,z) = (1/2) [cosh(cqz)/cosh (cqd) + 

sinh (cqz)/sinh (cqd)]. (69) 

It can be verified that g B {x, z; x , z ) satisfies Eqs. (21), (23), and (26) as required. Finally, 
we derive the boundary function by taking the PFT of Eqs. (28) and using Eq. (34) as follows: 

O') = V2/v Eq^ x P0Qx)<P h (q,z;x B ) 9 (70) 

where 

c P h (q,z;x B ) - g h (qf,z) V B (q), (71) 

and 

V B (q)= V 0 exp(-iqx B ). (72) 

We can verify by direct substitution that d> B (r) is the desired solution. 


3. BMGF for a composite solid 

Now we consider a 2D composite solid as shown in Fig. 1, consisting of the host H and an 
inclusion I. Our formulation is valid for an inclusion of arbitrary shape, but, for the sake of 
definiteness in this paper, we assume the inclusion to be an ellipse, which is a reasonably general 
2D shape. 

We calculate the BMGF for a single point source at r\i in H and induced sources in I. The 
induced sources are virtual sources introduced for the purpose of satisfying the boundary 
conditions at the interface as explained below. Since there is an inclusion, with properties 
different than the host, there must be boundary conditions to be satisfied at the interface of H and 
I. We consider the limiting case when the conductivity is infinite for a metallic inclusion. In this 
approximation, the temperature for thermal problems or the potential for electrostatic problems 
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must be constant over the entire interface. Without loss of generality, we assume this constant to 
be zero. Thus, the boundary condition at the interface is 

Hr,) = <t> 0 (r,) + <t> B (r,) = 0, (73) 

where Vi = (x 7 ,z 7 ) is a point at the interface and I is a running index. The function d> B (r 7 ) is 
known as given by Eq. (70). We choose a total of q discrete points over the interface, labelled 
1=1,2, — q. at which Eq. (73) is prescribed. 


Ideally, Eq. (73) should be satisfied at all points on the interface. In practice, we satisfy this 
equation for q points and choose q to be large. Note that Eq. (73) is for the total solution and not 
just for the GF. In writing this equation, we have used the same form of the total solution as 
given in Eq. (10). 

By using Eqs. (33) and (73), we obtain the following interface condition for the GF: 

G(r 7 ;r M ) = O 0 (r 7 ) = - <£ B (r 7 ), (74) 

where, by using Eq. (70), 

^bO/) = 1 / 2 m Iqexp(tqx 7 )d) 7l (q,z 7 ;x B ). (75) 

In order to satisfy the boundary/interfacial condition given by Eq. (74), we need another 
independent homogeneous solution. The GF method provides a useful technique for generating 
independent homogeneous solutions by introducing virtual sources or charges [1]. It is similar to 
the method of images [11, 25] in electrostatics or virtual forces [41, 43] in elastostatics. In this 
method, we introduce virtual sources just outside the boundary of the solution domain. If the 
virtual sources are denoted by the unknown function f (r), 

f(r) = fj 8(r - r m ), (76) 

where 

l r m| = M - 8 , 
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(77) 



I'm ^ r i j 


( 78 ) 


and 8 — 0 + in the limit (see Fig. 1). The positive sign of 8 ensures that the locations of the virtual 
sources are just outside the solution domain H. It also ensures the inequality (78) holds, which 
avoids the singularity in the GF calculation. For the purpose of calculating the limits, it is 
understood that (i) the point r = IT is exactly at the interfacial boundary, (ii) it is approached 
from inside the solution domain for calculating the boundary value of the solution, and (iii) the 
virtual sources are applied on the other side of the interfacial boundary, as prescribed in Eq. (76). 

This procedure ensures that the solution given by the virtual sources is the homogeneous 
solution, as described in the previous section after Eq. (41). It also determines the sign of (z-Zj), 
which has an important effect on the value of the GF and its derivative. 

The q arbitrary parameters fi are the unknown virtual source parameters, each uniquely 
associated with the corresponding point rj, according to Eq.(76) . The arbitrary parameters fi are 
determined by imposing Eq. (73) for the metallic inclusion. The solution corresponding to the 
virtual sources is then given by Eq. (27) with p(r r ) replaced by fb at r / = r H [. 

Including the contribution of the virtual sources in Eq. (54), the solution in domain H for a unit 
source located at r = r M becomes 


G(r; r M ) = G p (r ; r M ) + G p (r; r HI ) f, + G h (r ). (79) 

Note that the first tenn in Eq. (79) is to be multiplied by Fm as in Eq. (32). As in the previous 
section, though G h (r) does not explicitly depend upon r jyj and Vj , the arbitrary constants in 
the homogeneous solution will depend upon these parameters through the boundary conditions. 

To write the GF in terms of the PFTs in our semi-discrete model, we recall the following 
definition from the previous section: 


G(r;r M )= 1 / 2 ^Iq ex P (Lqx) T t (q,z; z M ) , (80) 
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where Y t (jq,Z) z M ) is the PFT of G (r; r M ). As in the previous section, the first, second, 
and the third terms on the RHS of Eq. (79) are calculated by using the following formulae: 


GpO; r M ) 

G P (X) r HI ) 

and 


^h( r ) - 



SqexpftqCx -X M )] T p (q,z-z M ), 



I^exp [iq(x - x HI )] T p (q,z - z HI ), 


exp(tqx) f h (q,z), 


(81) 


(82) 


(83) 


where ( q , z ) is given by Eq. (62) and ( q, z) is given by Eq. (63) with different arbitrary 

constants. 


We detennine the arbitrary constants in the homogeneous solution T \ (q, z) by applying Eq. 
(24) and then determine the virtual sources fj by using Eq. (74) for the interfacial conditions. 
For notational convenience, henceforth, we will denote the subscript HI on x and z in Eq. (82) 
simply as I. It should not cause any confusion about the location of the virtual sources relative to 
ij, which is defined by Eq. (77). 

The PFT of Eq. (79), which is a function of q and z, is given by 


r t (q,z; r M ) = exp ( - tqx M ) T p (q,z - z M ) + 

I 7 exp (-tqx 7 ) F p (q,z - z^ ft +T h ( q,z ). (84) 

Using Eqs. (8) and (9) along with Eqs. (47) and (48), we require 

r t (q,d;r M ) = 0, (85) 


25 



and 


F t (q,-d;r M ) = 0. 


( 86 ) 


Using Eqs. (62), (63), (83) - (86) yields the following expression for 


r\(q, z) = [e 1 (q)cosh(cqz)/cosh(cqd) + 

e 2 (q) sinh (cqz)/sinh (cqd)], (87) 

where 

e i(q) = - exp (-t qx M ) T sM (q) 

- £/ exp(—t qx 7 ) r s7 (q) / 7 , (88) 

e 2 (q) = - exp(-tqx M )r mM (q) 

— exp(—t qxj ) Tm/Cq) //, (89) 

r 5 /(q) = [r p (q, d - z 7 ) + r p (q, -d - z 7 )]/2, (90) 

r m /(q) = [r p (q,d - z 7 ) - r p (q, -d - z 7 )]/2 , (91) 


and r s7W (q) and T c7M (q) are given, respectively, by Eqs. (66) and (67). 

It can be verified that Y t in Eq. (84) with ei(q) and e 2 (q) given by Eqs. (88) and (89), will satisfy 
Eqs. (85) and (86) for all values of f 7 . Note that with as given by Eq. (87), fy and hence G 
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will satisfy Eq. (24) for all values of r M , r It and / 7 . The homogeneous solution depends upon 
the virtual sources fj that are yet undetermined. They will be detennined by imposing the 
interface condition given by Eq. (73). 

We write fy(c), z) in Eq. (87) as a sum of known and unknown terms as follows: 
r h 0 q,z ) - exp(—t qx M ) T hM (q,z) + 

I/exp(-t qxj ) r w (q,z)fy , (92) 

where 

r hM {q,z) = - [ r sM (q)cosh(cqz)/cosh (cqd) + 

r mM (q) sinh(cqz)/sinh (cqd)], (93) 

and 

T hI (q,z)= — [ r s7 (q) cosh(cqz)/cosh (cqd) + 

rm/ (q) sinh(cqz)/sinh (cqd)]. (94) 

We consolidate the above results into the following form that gives the GF at all points r in the 
solution space: 

G (r; r m) = (r; r M ) + £ 7 W x (r; r 7 ) / 7 , (95) 

where 


W M (r; r M ) = G p (r ; r M ) + ^V 2iV J Iq exp[tq(x - x M )] r^fy^z), (96) 

Wi(r; r 7 ) = G p (r; r 7 ) + fV 2 )v) ex PN(* - x 7 )] r w (q,z), (97) 


where we have used Eqs. (83) and (92). The subscript I on W is not a running index. It is used to 
indicate that this term is a characteristic of the interface in contrast to the subscript M which 
denotes the Mindlin source. 
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Once again, it is instructive to verify that G(t", r M ) in Eq. (95) satisfies Eqs. (17) and (22) for all 
values of r M , r h and fj. Note that the second term on the RHS of Eq. (96) and both the tenns 
on the RHS of Eq. (97) are solutions of Eq. (31), the homogeneous equation, in the solution 
space, because, in view of Eq. (77), r t is outside the solution space. 

To determine f,, we use the interface conditions given by Eq. (74) for G(r; r M ) at the r| values of 
r =1*! . This gives the following Ni equations for J =1,2 ... Ni, in Ni unknown parameters fj 
(1=1,2 ..Nj). 


- <l> B (r y ) = M/„(r y ; r„) + W,{r f r,) /, . (98) 

We can now write Eq. (98) in Ni -dimensional vector space of the coordinates of the virtual 
sources (b = 1,2, — Ni) as follows: 

Qj = 'Z,Vjjfi (99) 

or, in matrix notation, 

Q=vf (loo) 

where Q and / are vectors (column matrices) and V is a square matrix. Their matrix elements 
are given below: 

Qj = - 4> B (r y )-W M (r y ;r M ) (101) 

V,J= Wi(r ; ;r,) (102) 

Finally, we obtain the following result for the virtual source parameters: 

f=V~ 1 Q. (103) 


Using Eq. (103) for f in Eq. (84), we obtain the final GF for all points in the region H for a 
composite containing a metallic inclusion. 
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4. Application to phosphorene nanocomposite containing a metallic inclusion 

As an illustration of the semi-discrete GF, we apply it to a nanocomposite consisting of black 
phosphorous or phosphorene (matrix, H) and a metallic inclusion (I). We assume that the two 
materials have a perfectly planar lattice structure with smooth boundaries. Hence their thermal 
or electric response can be represented by the 2D Laplace/Poisson equation. 

Actually, phosphorene has a puckered layered structure and is not exactly planar like graphene. 
Strictly speaking, a mathematical model for nanomaterials should account for their discrete 
atomistic structure. However, for modeling thermoelectric properties of devices for industrial 
applications, finer effects at the atomistic scales are not always needed. Even for rigorous 
atomistic modeling, it is useful to first make a preliminary estimate of the required design and 
characterization parameters using the continuum model. A rigorous atomistic model is 
computationally expensive and has some inherent physical as well as mathematical uncertainties. 
On the other hand, the continuum model has been used for more than two centuries. It has the 
important advantage that it can provide a reliable preliminary estimate, which can provide a 
starting point for more detailed atomistic calculations. This explains the strong interest in 
computationally efficient techniques for solution of the Laplace/Poisson equation even for 
nanocomposites. 

The analysis given in this section is valid for any 2D material or nanocomposite, subject to the 
assumption of 2D planar structure and the continuum approximation. In what follows, we 
assume, without loss of generality, that the off-diagonal element of the symmetric conductivity 
tensor C is zero. If not, the coordinate axes can be transfonned to make C diagonal. 

For phosphorene, we assume that C 33 = 110 W/m-K and C 11 = 36 W/m-K. These values are 
actually theoretical estimates [9] for the armchair and zigzag directions in black phosphorene. 
These two directions may not be along the two Cartesian axes in a real phosphorene crystal. 
However, in this paper, to illustrate the effect of anisotropy we assume these two values of 
conductivity along the X and the Z directions. Consistent with the definition of GF, being the 
response to a unit probe, we assume F M = 1. 
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The numerical results are presented and discussed in the next section. A summary of the 
numerical values of the parameters and definitions used in the calculations is given below with 
reference to Fig. 1. 

Host: 

P = + ct and /? = - Cl (104) 

C = Cii/C 33 = 0.3 (105) 

L=40, d=40 

Coordinates of the Mindlin source point, where F M is applied: r M = (0,0) for Fig. 2 and (-30,-20) 
for Fig. 3. 

Applied source strength at r M : F M = 1 
Dirichlet boundary condition at outer boundaries: 

G(0,d) = V 0 = 1. 

G(x, z) =0 at all other points (all values of x except 0) for z = +d and -d. 

Periodicity in the X direction: 

G(x+2nL, z) = G(x, z); 

where n = any integer (positive, negative, or 0), subject to the limits given in Eq. (44). 

Number of discrete points in the X direction between 0 and L (three sets): N= 100, 200, 400. 
z is a continuous variable. 

Inclusion: 

a = OA = 1; a = OQ = 0.5. This corresponds to eccentricity of the ellipse V[l-(a’/a) 2 ] = V 3/2. 

G at all points on the interface as given by Eq. (74). 

Number of discrete points at the interface where boundary value is specified: Ni = N*a/L 

5. Results and discussion of the convergence of our semi-discrete model 
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We first calculate the GF for homogeneous phosphorene without any inclusion. We specify the 
boundary conditions given by Eqs. (8) - (6). In order to check the numerical convergence of the 
GF, we calculate it for three different values of N- 100, 200, and 400. Values of GF along the 
perimeter of the ellipse in Fig 1 are shown in Fig. 2 for these values of N. A point on the ellipse 
is labelled by the angle 0 of its radius vector from the positive X-axis. Thus 0 = 0 (and 2n) is the 
point S in Fig. 1. Similarly, 0 = rc/2, n, and 3n/2 are, respectively, the points Q, P, and T in Fig. 1. 
In Fig. 2, the GF values for N=200 are shown by the solid curve. The values for N =100 and N = 
400 are shown as discrete points on the graph. We have done the calculations for many more 
points but much fewer points are shown in Fig. 2 for visual clarity. For perfect convergence, the 
GF should be, obviously, independent of N and the points for N=100 and 400 would he on the 
curve. In Fig. 2, the points for NQ=400 lie almost exactly on the curve for NQ=200, which 
shows excellent convergence at N=200. Even points for N=100 he quite close to the curve. On 
the scale of the graph in Fig. 1, they are not visually distinguishable from the curve N=200. We 
infer from Fig 2 that the GF converges very well for N= 200 and quite well even for N=100. 

The convergence test in Fig 2 is more stringent than it might appear. Usually the convergence is 
tested by calculating the values of functions at the same space points using a variable number of 
discretization points. This test does not guarantee that the values of the function at intennediate 
points have also converged. In contrast, the GF in Fig. 2 for different N has been calculated at 
different space points but of course on the perimeter of the same ellipse shown in Fig. 1. The fact 
they all lie on the same GF curve clearly shows that the values of the GF have satisfactorily 
converged at all points (original mesh points as well as intennediate points) on the ellipse. 
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Angle from X axis 


Fig. 2: Convergence of GF with respect to number of discrete elements in the X 
direction for a homogeneous solid without inclusion. Figure 2 shows the value of the 
GF at different points on the ellipse in Fig. 1, identified by the angle of the radius 
vector from the positive X axis (see text). The ordinates Q, P, and T are the same as 
labelled in Fig 1 and correspond, respectively to the angles tc/2, n, and 3n/2. The 
source is at (0,0). The solid curve is for NQ= 200. Values for NQ=T00 and NQ = 400 
are also shown. 
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Using the GF for the homogeneous solid as just calculated, and the parameters as given in the 
previous section, we calculate the BMGF for the composite shown in Fig. 1, with the boundary 
and the interface conditions given by Eqs. (8), (9) and (73) for a metallic inclusion. The 
calculated values of the BMGF are shown in Fig. 3. For points inside the inclusion, the value of 
the GF is taken to be equal to that at the interface, which is 0. The calculated BMGF can be used 
to calculate the full solution for any integrable distribution of Mindlin sources in the host region 
and Dirichlet values at the outer boundaries. 
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1 0.000 
0.2150 
0.4300 
0.6450 
0.8600 
1.075 

1 1.290 
1.505 
1.720 
1.935 
2.150 


Fig. 3: Calculated BMGF shown as a function of x and z for phosphorene containing a 
metallic island at the center (see Fig. 1). The big peak is at the source point (Mindlin point) 
and the small peak at the outer boundary z=40 is where the Dirichlet condition given by Eqs. 
(11) and (34) (Boussinesq point) is specified. Compared to the height of the peaks, the BMGF 
at other points is almost zero. The values inside the elliptic inclusion are constant, which is set 
to zero in the calculations. Flowever, in the figure, the constant is given a negative value in 
order to show a better contrast with the values just outside the inclusion. 


Figure 3 provides the desired map of the response of the solid over all the space points. We note 
that, as expected, the BMGF has a strong singularity at the Mindlin source point and a relatively 
weak singularity at the Boussinesq boundary point where the Dirichlet condition has been 
imposed. Moreover, the shape and size of the inclusion are also apparent in Fig. 3. This shows 
that, at least in principle, the measured GF can be used for a partial characterization of the 
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composite. Further, the analytic nature of our model also indicates the feasibility of solving the 
inverse problem of determining GF from the measured response of a composite. 

In our semi-discrete model, the solid is defined only at discrete values of x given by Eq.(43). As 
mentioned before, this feature is similar to the Bom-von Kannan model of a crystal lattice except 
for the following important difference. In the Bom-von Kannan model, the space coordinates of 
the atoms are fixed but the actual dimensions of the lattice change with different values of N. In 
our model, the dimensions of the solid are fixed but the positions of the elements are different for 
different values of N. If the objective is to calculate the GF at a specific point in space for 
different values of N, then the X-coordinate of that point has to be scaled appropriately. 

We can estimate the convergence efficacy of our method for solving the Laplace/Poisson 
equation for a homogeneous solid as follows: A numerical technique essentially approximates a 
continuum space by discretizing it, as for example in the boundary element and the finite element 
methods, and in numerical integration techniques. This approximation breaks the original space 
into a set of discrete elements, with the solution over each element approximated by a convenient 
interpolating function. There are two main sources of error in such techniques: (i) the shape and 
size of the elements may be such that they don’t fill the whole space exactly, and (ii) the error 
due to the approximating function not being an exact solution. 

In our model, the error (i) is not present because the elements are continuous in the Z-direction. 
The shape of each element can be simulated exactly by specifying the functional dependence of z 
as z = fix). So there are no spaces left because of mis-fitting elements. Regarding the error (ii), in 
general, the error or the uncertainty due to an approximating function in discrete elements can be 
assumed to be depend inversely (but not necessarily linearly) on the number of elements or 
directly on the area of each differential element AAh, where Ah is the area of the solution space. 
We can write 

Ah = L 1 L 3 (106) 

where Li and L 3 are the relevant lengths in the X and the Z directions. In the present case of a 
square solid, Li = L 3 . The full discretization consists of discretizing the solution space in both the 
X and Z directions. Consistent with this assumption, the error or the uncertainty due to one¬ 
dimensional discretization in the i th direction (i=l or 3) can be assumed to be equal to the length 
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of the differential element AL; in that direction. If N; denotes the number of elements in the i- 
direction, 

AL, Lj/Ni. (107) 

Using Eq. (106), the relative error 8 d in the solution due to the full discretization of space can 
then be estimated as follows: 

£ D AA h /A h = ( Li AL 3 + L 3 ALi)/LiL 3 = ALi/Lj + AL 3 /L 3 . (108) 

Eq. (108) does not really give a quantitative estimate of the error because the actual error will 
depend upon the nature of the function which is being approximated. At best, it gives a 
parameter that can serve as a qualitative measure of the error. In a comparison between two 
systems, the numerical error in a system with a larger value of 8 d is likely to be more than a 
system with a lower value of 8 d- We will refer to 8 d as the relative error only with that 
understanding. 

Now we estimate the relative error in our semi-discrete model. In this model the solution space is 
discretized only along Li. The solution in z is obtained exactly, which does not depend upon the 
element size. Hence AL 3 = 0 in Eq. (108). Thus, 8s, the estimated relative error in our semi¬ 
discrete model is given by 

8s AA h /Ah = ALi /Lj. (109) 

Assume that the size of the elements in a fully discrete calculation approximately same in the x 
and the z-directions or Thus, we see from Eqs. (108) and (109) that the computational error 
introduced due to discretization in the semi-discrete model is of the same order of magnitude as 
in the fully discrete model. At first glance, this inference may appear to be surprising, but there is 
an important difference. 

To keep the errors in the two models comparable, that is, if we require, 

£d £s (HO) 


then, from Eq. (107), 


36 



N, N 3 N 


( 111 ) 


2 

We see from Eq. (Ill) that a fully discrete model would require computations at N 1 N 3 N" 
elements for the error to be comparable to our semi-discrete model. In contrast, our semi-discrete 
model would require computation only at N points. Since N 10" - 10 or even larger, the 
computational efficiency of our approach in terms of the computational cost should be one or 
two orders of magnitude better than a conventional fully discretized approach. 

6. Conclusions and general discussion 

In this paper, we have described a semi-discrete model for calculating the GF for the 2D 
Laplace/Poisson equation for the modeling of the thermal/electrical response of 2D 
nanocomposites. The method is computationally very efficient and accounts for anisotropic 
conductivity, which may be important for modeling 2D materials. 

The main features of our semi-discrete model and the calculation of GF using PFT are 
summarized below: 

i. Out of the two field variables (x and z), only one variable (x) is discretized and is restricted to 
prescribed discrete values. The other variable (z) is continuous. In 3D, two variables, x and y, 
will be discrete and z will still be continuous. The discrete variable is allowed only certain values 
in a finite range, in accordance with the periodic boundary conditions. This is similar to the 
Born- von Kannan model of lattice dynamics, in which the conjugate Fourier vector is also finite 
and discrete and is confined to the first Brillouin zone. This makes it easy to impose periodic 
boundary conditions that are automatically obeyed. The integral over the discrete space and 
conjugate variables are expressed in terms of discrete sums with well-defined orthonormal 
properties and convergence characteristics. 

ii. The continuous nature of z makes it possible to reproduce all inclusion shapes exactly. If all the 
variables have to be discretized as in the conventional numerical models, the model would 
necessarily leave out some space between discrete points. 

iii. As shown in Sec. 2 and Sec 3, the use of the GF allows modeling to be done in modules of 
increasing complexity. This avoids wasted computational effort since the result of any module 
can be used in the next module. For example, for the problem in this paper, we first calculated 


37 



GF for a homogeneous solid with only periodic boundary conditions, which was the 0 th module. 
The GF for this module is given by G p in Eq. (49). In the next module, we introduced boundary 
conditions in the Z-direction. In the final module, we simulated a composite solid with 
prescribed conditions at the interface. The GF for each module obeyed the boundary conditions 
of the previous module and only the conditions prescribed for the next module had to be 
imposed. 

iv. Our model is semi-analytical, in the sense that the Fourier integral over one of the variables is 
done analytically and, therefore, exactly. Because of this fact the convergence efficiency of our 
technique for a homogeneous solid should be at least an order of magnitude better than the 
techniques in which both the variables have to be discretized. For nanocomposites, the 
computational efficiency of our technique is discussed in the next paragraph. 

v. To estimate the relative computational advantage of our technique for a nanocomposite, we note 
that there are two main steps in the calculation: one, the calculation of the homogeneous BMGF 
and two, the calculation of the BMGF for the nanocomposite. We have shown in the previous 
section that the numerical convergence in our semi-discrete model can be an order of magnitude 
faster than the conventional fully discretized models. For modeling a nanocomposite in step 2, 
we need to invert a matrix for calculating the virtual sources in Eq. (103). The order of the matrix 
to be inverted is Ni x Ni, where Ni is the number of points on the interface at which the interface 
conditions are to be satisfied. Thus, for calculation of the BMGF in our model, we need a 
numerical mesh only at the interface between the inclusion and the host. The number of points 
on this mesh is obviously much smaller than in the conventional finite element or boundary 
element method in which a numerical mesh is needed either on the whole solid or at all the 
boundary points. Since the computational effort in matrix inversion is O (N ) (for an N x N 
matrix), reduction in the number of mesh points at which a matrix needs to be inverted results 
into a major economy in the computational budget. 

We can estimate the computational advantage in the present case as follows. In case of the BEM, 
the number of mesh points is proportional to the total length of the boundaries. The outer 
perimeter of the solid H is 2(L+d) (see Fig. 1). If the perimeter of the ellipse is 2e, the total 
length of the boundary in the BEM on which the numerical mesh needs to be created is 
2(L+d+e). In contrast, the corresponding length in our model is just 2e. For e = 1/10 of the outer 
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boundary, the order of the matrix to be inverted in our model is approximately 1/10 of that in the 
conventional BEM. Since the matrix inversion is an 0(N ) process, a crude estimate of the 
saving in CPU effort can be as much as by a factor of 10\ 

In case of finite element modeling, the number of mesh points is proportional to the surface area 
of the solid. This factor is much more than even the BEM. However, it is not a valid estimate of 
the computational advantage, because the relevant matrices in the FEM are sparse and their 
inversion requires much less than the 0(N 3 ) process. A better estimate will be that the matrix 
inversion estimate in the FEM is about the same as the BEM. Again, FEM has been a very 
popular technique that has benefitted by excellent refinement in the field. Moreover, FEM has 
some additional advantages such as including the non-linear terms whereas the GF method is 
essentially a linear method. 

It is admittedly difficult to put a precise number for the numerical advantage of our technique for 
modeling composite materials. However, in view of the above discussion, we can claim at least 
qualitatively that the proposed technique should result into a substantial economy in the 
computational effort - even by one or more order of magnitude. It should result into a robust and 
reliable algorithm for modeling of nanocomposites. The actual economy will obviously depend 
upon the specific problem. An attractive possibility may be to combine the GF method with the 
conventional BEM or FEM methods. 

Finally, we reiterate that the GF gives the response to a probe, and is, therefore, a measurable 
physical characteristic of a composite nanomaterial [1]. A reliable and robust calculation of the 
GF, as presented in this paper, should be useful in designing and interpreting an SPM based 
characterization technique for 2D materials. 


Acknowledgements 

The authors thank Drs D.T. Read and A. Smolyanitsky for their valuable comments on the 
manuscript. 


39 



Appendix A: Evaluation of integrals for GF and the flux 


In this appendix, we evaluate the singular integral on the RHS of Eq. (35) , which is defined in 
the Cauchy sense. For notational brevity in this Appendix, we assume that the source is located 
at r=0. For a source at r M , replace r by r- r M . 

We introduce a small imaginary part in the denominator of the integrand in Eq. (35) as follows 

[40] : 

G(r) = (1/2 tt) 2 J_ + ” lim £ ^ 0 [H*(k) - ie]' 1 exp[tk.r] dk. (112) 

From the real part of the inverse on the RHS of Eq. (112), we obtain 

G P (X) = (1/2tt) 2 P/_ + “ [W(k)]- 1 exp[tk.r] dk, (113) 

where the subscript p denotes the particular GF introduced in Eq.(30), P denotes the Cauchy 
principal value, and we have used the following relation, which is valid for any real scalar 
variable 

lim £ ^ 0 ^ = pQ) + tit6(5). (114) 

In our semi discrete model, as given in Eq. (45), k = (q, k 3 ) and r = (x, z). Using Eqs. (40) and 

(41) , we write Eq. (113) in the form 

G p (r) = (1/2Af ) SqT p (q,z) expftqx], (115) 

where we have replaced the integral over ki by a summation over q in accordance with Eq. (46), 

Y v (q,z) = (^)P/_ + ”G p (q,/c 3 ) exp(t/c 3 z) d P 3 , (116) 

and 

6 p (q,fc 3 ) = ['P(q,/c 3 )]- 1 . (117) 

From Eq. (40) 

<P(q, fc 3 ) = C 33 (k 3 -q/S)(k 3 -qfi)]. (118) 
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In Eq. (118), q/3 and its complex conjugate q/3 are the two roots of the following equation for 

k 3 , 

fcf + 2(C 13 /C 33 )qk 3 + (C n /C 33 ) q 2 = 0. (119) 

We denote the real and imaginary parts of P, by p R and Pi respectively, 

P=P R +i0 b ( 120 ) 


with the assumption 

Pi >0. (121) 

The two roots q/3 and q(3 correspond to zeros of the function ^(k) in Eq (40). The two roots 
of Eq. (119) must be either both real or complex conjugate of each other. This is because the 
LHS of Eq. (119) is real. In general, the off-diagonal element of the conductivity tensor C is 
smaller than the diagonal elements so that Cn < C 11 C 33 . Hence the two roots are complex and 
must be conjugate of each other. Without loss of generality, we choose P such that Im (P)>0. 

First, we consider the integral in Eq. (116) for q > 0. The result for q < 0, can be obtained in an 
analogous way. Equation (117) can be written in the following alternative forms: 

GpCq.fc,) = 1/(2 iC 339 /? ; ) [l/(k 3 - qp)- l/(k 3 - q/S)], (122) 

or 

GpCq.fc,) = l/(2iC 33 p, k 3 ) [p/(k,-qp)-p/(k 3 -qp)]. (123) 

To carry out the integral in Eq. (116), we use the Cauchy’s residue theorem [43]. We use the 
form given in Eq. (122) and choose an anticlockwise semicircular contour. For z>0, we choose 
the contour in the UHP (upper half plane), as shown in Fig 2. The only pole that will contribute 
to the Cauchy integral is the one in the UHP, which, in view of Eq. (122), is at k 3 = q(3. This 
gives the following result: 

r p 0z,z) = exp(t^z) H(z)/(2C 33 ^ / ), forz>0 (124) 

where H(z) is the Heaviside step function, which is 0 for z<0 and 1 for z>0. 
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For z <0, we choose a clockwise semicircular contour in the LHP. In this case only the pole in 
the LHP would contribute to the integral. It gives the result 

T p (q,z) = exp(i^qfz)H(-z)/(2C 33 ^ / ) forz<0. (125) 

In deriving Eq. (125), we have followed the standard sign convention of anticlockwise contour 
being positive. Eqs. (124) and (125) can be written in the following concise fonn: 

Y p (q,z)= exp(ip R qz- faq\z\) /(2C 33 qfa). (126) 

The limiting cases of z = 0, q = 0 have to be treated directly and separately keeping in mind that 
the result will depend upon the direction in which the limit is approached. It can be verified by 
direct substitution that G p (r, r M ) given by Eq. (115) satisfies Eq. (30) and <l> 0 (r) given by Eq. 
(27) satisfies Eq. (1). 

The imaginary part of Eq.(l 14), gives the homogeneous GF [40] G h , which satisfies Eq. (31) 
for any r M- Using Eq (114) in Eq. (112), we obtain from the imaginary part 


G h (r ) 6[¥(k)] exp[tk.(r)]. (127) 


The constant of proportionality in Eq.(127) is arbitrary and is detennined by imposing the 
boundary conditions. 


Note that in Eq. (127), it has been assumed that the Mindlin source is located at =0. It can 
be verified by direct substitution [40, 41, 43] that G^(r; r^j ), as given by Eq. (127) for r 
replaced by r — r M is a homogeneous solution that satisfies Eqs. (21) or (31) for any V M . 
For this purpose, we use the following relation 

?6(0 = 0, (128) 
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valid for any real variable The constant of proportionality in Eq. (127) is arbitrary and is 
determined by imposing the boundary conditions. It is an interesting feature of the GF method 
that it can give both particular and homogeneous solutions in a single representation. 

We also quote the following identity for the delta function of any real function F in Eq. (127): 

5[F(?)] = S( ? - W/F'Ro ), (129) 

where 

F '© = dF/d?, (130) 

and is a solution of 

F(5o) = 0. (131) 

Note that F is real but can be complex. We infer from Eqs. (127) and (129) that G^ will be 
nonvanishing only for those values of k for which ^(k) = 0. 

Using the above properties of the delta function, we write the homogeneous GF in the same form 
as Eq. (115) for our semi-discrete model using the relations given in Eq. (45). As in Eq. (115), 
we write: 

Gft(r) = (l/2W) I, q T' h (q,z) exp(tqx), (132) 

where, by using Eq. (131), 

r h(q>z) = I,k s3 6s (q) 5[$(q,/c s3 )] exp(t/c s3 z), (133) 

and /c 3 = k s3 is a solution of ^(q, /e 3 ) = 0 or Eq. (119), and 0 s (q) is an arbitrary 
constant. However, it can be function of q and, in general, will be different for different k s3 . In 
view of Eq. (127), each term in the sum on the RHS of Eq. (133), and therefore the sum, will 
satisfy Eq. (31), and will be a homogeneous solution. In general, there will be two values of k s3 , 
which are the same as the two roots of Eq. (119). 
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Fig. 4. The semicircular contour in the upper half plane (UHP) for 
Cauchy integration in Eq. (116) for q>0 and z >0. The two poles are 
shown as thick dots. In this case, only the pole in the UFIP contributes 
to the integral. In this figure /3 R ^ 0. 
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